Theory of correlations between ultra-cold bosons released from an optical lattice 
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In this paper we develop a theoretical description of the correlations between ultra-cold bosons after free 
expansion from confinement in an optical lattice. We consider the system evolution during expansion and give 
criteria for a far field regime. We develop expressions for first and second order two-point correlations based 
on a variety of commonly used approximations to the many-body state of the system including Bogoliubov, 
meanfield decoupling, and particle-hole perturbative solution about the perfect Mott-insulator state. Using these 
approaches we examine the effects of quantum depletion and pairing on the system correlations. Comparison 
with the directly calculated correlation functions is used to justify a Gaussian form of our theory from which 
we develop a general three-dimensional formalism for inhomogeneous lattice systems suitable for numerical 
calculations of realistic experimental regimes. 
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I. INTRODUCTION 



In the past decade noise correlation analysis, analogous to the pho- 
ton correlations observed in the landmark experiments of Hanbury- 
Brown and Twiss fj]], has been applied to ultra-cold atom experi- 
ments. In atomic systems such measurements can be used to reveal 
information about interaction-induced (i.e. many-body) correlations 
between the atoms. In addition, the dramatic differences between 
Bose and Fermi statistics have been clearly demonstrated with neu- 
tral atoms (e.g. see U). 

A wide range of atomic physics experiments examining various 
aspects of correlations have been conducted. The first experiments 
by Yasuda et al. (J observed atom bunching using a neutral beam 
of (bosonic) neon atoms. In quantum degenerate Bose gases local 
third-order correlations have been inferred by measuring three-body 
decay rates [J, |5|], and first order coherence has been studied using 
matter wave interference @, 0] an d Bragg spectroscopy 0|U. Of 
most interest for the investigation in this paper has been the recent 
experimental progress in the spatially resolved measurement of sec- 
ond order correlations in both bosonic and fermionic ultra-cold gases 
(1 0, 0- Two general approaches are used to make 
these measurements. One approach involves directly counting atoms 
(H 0, [Till, while the other uses absorption imaging to measure the 
density ^13 , , , 0, 01 . The applications of these measurements 
have included Bose and Fermi gases in harmonic traps ffllicWlHliUl. 
and in optical lattices 001151 0. 

In this work we are concerned with the theoretical formalism for 
the spatial noise correlations of an ultra-cold Bose gas after expan- 
sion from an optical lattice, relevant to the experiments reported in 
Refs. I"r2l0. [l6ll . Initial theoretical work on this subject was pro- 
vided by Altman et al. fnll . who predicted the noise correlations 
using a perfect Mott-insulator approximation (i.e. neglecting tunnel- 
ing) and assuming a simplified form for the single particle expan- 
sion. Subsequent experiments in the Mott-insulator regime fl2l [T^l 
verified those predictions, in particular the periodic bunching peaks 
in the noise correlations. Several recent theoretical proposals have 
built on that framework and investigated the use of noise correla- 
tions in characterizing many-body states produced in optical lattices 

SO HO, EE [Hdl. This line of research provides an interesting 
new avenue for investigating the effects of interactions which com- 
pliments the other techniques available such as direct density imag- 
ing 000L Bragg HE 0,00 00 and Raman spec- 
troscopy 1313, . 

The basic organisation of the paper is as follows. In the remainder 



of this section we introduce the system of interest and give an intro- 
duction to how the far-field correlations are determined. In Sec.lTTlwe 
derive the properties of single particle expansion from the lattice and 
use this to derive a simplified far field form and its validity condi- 
tions. We discuss the correlation function formalism and its relation- 
ship to the Bose-Hubbard Hamiltonian in Sec. [Ill] The main results 
for expanded correlations functions of a ID lattice system are devel- 
oped in Sec. llVl using a variety of theoretical approaches. In Sec. [VI 
we introduce a Gaussian approach which we justify by comparison 
to the earlier results. In Sec. lVII we extend this Gaussian approach to 
a general 3D theory for the inhomogeneous lattice system, and then 
conclude. 



A. Optical lattice 



Consider a system of bosonic atoms in an optical lattice, described 
by the Hamiltonian 



H = J d 3 xft(x.) \J?L + J2 V sm 2 (k Xj ) + V oxt (x) 



+ ^t(x)^(x) 



(1) 



where ?/j(x) is a bosonic field operator, V cx t describes any external 
potential (typically harmonic) and Uo = 4,-Ka s h 2 /m characterises 
the binary interactions between the particles, with a s the s-wave scat- 
tering length. The lattice is taken to be simple cubic, with a = n/k 
and b = 2k the lengths of the direct and reciprocal lattice vectors 
along each direction, where k is the wavelength of light used to 
produce the lattice. This lattice is hence of separable form, of the 
type produced in experiments with three sets of orthogonal counter- 
propagating light beams, and the depth along each direction (i.e. Vb) 
is assumed to be the same. The theory we present here can be easily 
extended to more general lattice configurations, but we restrict our 
attention to this case for notational simplicity. We define the quan- 
tities lor = hk 2 /2m and Er = Fllor as the recoil frequency and 
energy respectively. 
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B. Measurement of expanded correlations 

Due to the external confinement (V ox t) the equilibrium state of 
this system will be of finite size Ma, where M is an integer and 
a is the lattice site spacing introduced above. Typically in experi- 
ments M < 100. It is of interest to probe the in situ correlations of 
this system, however imaging limitations make this difficult and it is 
necessary to turn off the lattice and external potentials and let the sys- 
tem expand first. Expansion itself is a rather complex process. The 
interactions between particles and high momentum components aris- 
ing from the tight lattice confinement, mean that spontaneous s-wave 
scattering into unoccupied modes will generate additional noise and 
correlations in the system. Here we neglect these collisional pro- 
cesses and assume that upon release the atoms freely evolve 16711 . 
After a short time of expansion the system density is sufficiently low 
that this is a good approximation, however this could be arranged 
experimentally by use of a Feshbach resonance. 

Within this ideal expansion approximation the procedure for de- 
termining the expanded noise correlations is as follows: The atoms 
are initially in the optical lattice and are described by some (inter- 
acting) manybody state |$). The confining potentials are suddenly 
removed at t = and the system freely evolves according to the 
propagator U(t) = e - iKt/h , with K = J d z x ft (x)[p 2 /2m]ip(x). 
The mean system density at time t is given by 

(n(x)) t ee (*\rf(t)ft(x)$(x)U(t)\*). (2) 

Experimentally the density can be measured using absorption imag- 
ing, however in any given measurement shot noise will be present 
and l(2j corresponds to the average of many such measurements. 
It has been demonstrated fl5l l24ll that under experimental condi- 
tions the atomic noise can be made to dominate the photon noise 
in the measurement process so that density correlations can be 
inferred from the noise analysis. In particular, the covariance 
C(x,x') = (n(x)ri(x')} — (n(x)} (n(x')} can be measured using 
absorption images (6§j], while the correlation function G' 2 ' (x, x') = 
(ip (x)-i/v (x')-i/j(x)-i/>(x')} is typically measured when direct atom 
counting is employed. 

It is of interest to know the relationship of the noise correlations 
signal to the in situ properties of the system. It is usually assumed 
that for large t we can make & far field approximation whereby the 
density expectation value becomes propotional to the momentum dis- 
tribution of the initial state, i.e. 

777 

(«( x )> ~ ^<«Q(x)}, 0) 

where 

Q(x) = mx/lit (4) 

relates the i>7 situ momentum hQ of the field to the final observation 
position x. Another important length scale we introduce at this point 
is the spatial separation of diffraction peaks which occur when the 
system is coherent. This length scale is set by the reciprocal lattice 
vector and is given by 

h(t) = hbt/m. (5) 

We also note the work of Kolovsky |36|| . who has performed nu- 
merical simulations of the second order correlation function for a 
small ID lattice (in the pure Bose-Einstein condensate and perfect 
Mott-insulator limits) and observes interesting changes between the 
near-field and far-field correlations. 
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Figure 1: Schematic of atomic far-field expansion from a lattice: 
Consider an atom initially [t = 0) located at lattice position Rj, 
and after expansion time t is found at far-field position x (measured 
relative to the lattice center indicated by position 0). The initial lat- 
tice sites are indicated as grey filled-circles. 



II. SINGLE PARTICLE EXPANSION FROM A LATTICE 

A. Semiclassical derivation 

Consider Fig. [T] An atom initially in the lattice at location Rj 
expands in time t to reach location x. We define the displacement 
vector for this particle during time-of-flight as dj=x — Rj. Observ- 
ing that the average velocity of the particle must be v = dj /t, we can 
use the de Broglie relations to assign it a wavevector Qj = mdj / ht, 
and respective frequency w(Qj) = ft|Qj| 2 /2m. Thus we conclude 
that the total phase change occurring to this particle during expansion 
would be 

A^=Q J -d J -a;(Q J )i=^^. (6) 

In the far-field approximation where ]x| ^> |Rj|, we have |dj| = 
|x - Rj| 2 « x 2 - 2x • Rj, so that 

A4>j « - Q • Rj , (7) 

where Acf> — w(Q)t and Q = mx/ht, as in Eq. ©■ 

We now turn to a quantum mechanical treatment of this expansion 
and show that we get the same spatially dependent phase factor (i.e. 

Q Rj) 



B. Wannier state expansion 

To construct an accurate quantum description it is necessary to 
consider the wavefunctions inside the lattice. Of most interest to us 
is the case where the atoms are in the ground band, and the Wannier 
states, wo (x) form a convenient localized basis. For sufficiently deep 
lattices a tight binding approximation is appropriate, whereby we ex- 
panded the lattice potential to second order in the position variable 
about each lattice site minimum to obtain a local harmonic oscillator 
description with frequency 

WLatt = 2y/V / ErUr. (8) 

For simplicity of argument, in this section we restrict our attention 
to the one-dimensional case of a lattice of M sites, with locations 

Rj=ja-\iM-\)a, j = 0, 1, . . . , M - 1, (9) 
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and take M to be odd for symmetry. Because the Wannier states are 
separable our results and conclusions apply immediately along each 
direction in the 3D case. 



The ID Wannier states of the ground band can be approximated 
as harmonic oscillator ground states 



expansion 



w (x — Rj) 



-0--R,-) 2 /2x 



(10) 



Figure 2: Far-field expansion of atoms from a lattice: A tight-binding 
Wannier state at lattice position Rj expands for time t. 



where xq = ^/h/muj^^tt is the oscillator length. 



After free evolution for time t the harmonic oscillator approximation to the Wannier state evolves to (see Fig. [2} 



wo{x — Rj,t) = 



exp 



2[W{t)f 
I 



{x-Rj) 2 ht 
6XP 1 1 2[W(t)} 2 mx\ >( 



(11) 




t [sec] 

Figure 3: Evolution of Wannier state width W(t) in units of the direct 
lattice vector a for typical optical lattice parameters and two lattice 
depths: Vo = 5Er (black curve) and Vo = 50Er (grey curve). 
Results for 87 Rb in a lattice made from counter-propagating 850 nm 
light fields. 



where W(t) = x \J 1 + (ht/mx 2 ,) 
of the Gaussian envelope, and 9 = 
Refs. HH. 



is the time-dependent width 

i arctan(tfi/m:EQ). Also see 



1. Far-field limit 

For sufficiently long expansion times the Wannier states will ex- 
pand to be much larger than the total extent of the sites initially oc- 
cupied in the optical lattice. We can then approximate 



W{t) ~ ht/mxo, 



(12) 



i.e. the size of the system grows linearly in time proportional to the 
initial momentum spread of the Wannier state. We examine this ex- 



pansion for typical experimental parameters in Fig. [3] We see that 
after about 10 ms of expansion the individual Wannier wavepackets 
are of order 100 times the lattice site spacing. For the results we 
derive in this section, we assume that Eq. H2\ holds, but even more 
strictly that the system has expanded to be much larger than the spa- 
tial extent of the initially occupied lattice (i.e. max{i?j} ~ Ma). 
Thus we have the hierarchy of length scales 



(13) 



(14) 



W{t) > Ma > x Q . 

We also note that when Eq. i 1 2b is satisfied, we have 

W(t) 
h(t) ' 

so that the diffraction peaks will be further out than W(t) unless 
Vo > 16E R . 

The above length scales do not actually define the far-field limit. 
Instead, we define this limit as being when the following approximate 
form for the expanded Wannier states is accurate: 

w (x-R h t) ~ A{x,t)e- iQ(x)H \ (15) 
where the common (complex) amplitude of all Wannier states is 



A(x,t) 



7T 1/4 yPF(i) 



(16) 



and Q(x) is defined in Eq. 10). 

We can consider this simplification as two distinct approximations 
to the phase and amplitude of the state, and we now consider the 
nature of these approximations. In particular, in addition to the re- 
quirements d 1 3 fc we define the validity regime of the far-field regime 
in terms of the spatial region X (i.e. far-field positions satisfying 
\x\ <g X) over which the amplitude approximation is accurate, and 
expansion time T (i.e. t > T) for which the relevant phase errors 
will be small. The basic small parameter we will use is the ratio 
Rj /W which has a maximum value of Ma/W . 
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2. Amplitude approximation 

We first analyse the error in the amplitude of the Wannier state 
Eq. d 1 5 b associated with ignoring the initial offset in the origin of the 
Wannier state Jilt , i.e. in setting 



\w (x - Rj,t)\ — ► = \w {x,t)\. 



(17) 



To first order in Rj /W the relative error in making this approxima- 
tion is 



2W 2 (t)' 



(18) 



To maintain a small amplitude error we can restrict the spatial range 
over which Eq. d 1 5 b is applied, i.e. 



|x-| < X M {t) = 



2W{tf 
Ma 



(19) 



In practice this means that the the amplitude approximation is accu- 
rate for max{|s:|} ~ W. 



3. Phase approximation 

We now consider the error in the phase of the Wannier state 
Eq. d 1 5b associated with ignoring the quadratic dependence on Rj 
in the phase term of the full Wannier state dl lb . To first order in 
iij /W the relative error in this approximation is 



*?m = - Jh - 



2W(t)xo 



(20) 



Maintaining a small phase error between a pair of Wannier states 
separated by p-sites requires that the expansion time satisfies 



t » T p 



2 2 i 
p 7T 1 

4 LOR ' 



(21) 



Satisfying this condition for all Wannier states (i.e. taking p = M) is 
difficult for the parameters of current experiments, and may require 
going beyond the exp(-iQRj) approximation d 1 5b to analyze the 
expanded images. 

For the parameters used in Ref. lfl2ll where M ~ 80 we have 
that Tm ~ 0.8 s, whereas in the experiments an expansion time 
of 22 ms was used (also see Table |IJ. However, because they only 
examined the (incoherent) Mott-insulator regime, only short range 
correlations are important. For this situation the relevant quantity is 
Ti ~ 0.12ms, which gives the timescale for the approximate ex- 
pression d 1 5b to accurately predict the phase relationship between 
atoms expanding from neighboring lattice sites. In situations where 
long-range phase coherence is important it will be necessary to use 
much longer expansion times or go beyond d 1 5b in the analysis of the 
results. 




-0.5 0.5 

Q [units of b] 

Figure 4: The function F(Q) for M = 3 (black line), M 
(medium grey line) and M = 25 (light grey line). 



C. Quasi-momentum basis 

It is also useful to consider an approximate form for the time- 
dependent expansion of quasi-momentum basis states. These states 
are defined as a discrete Fourier transform of the Wannier states, i.e. 
as 



4> q (x,t) = — = ^2 e qR Jw (x - Rj,t). 

V 3=0 

Making the far-field approximation d 1 5b we obtain 

(f> q (x,t) w A(x,t)F M (Q{x) - q) , 
where we have defined 



HQ) 



sm(^MQa) 
\/Al sin(iQa) 



(22) 



(23) 



(24) 



We note that F(Q) is a peaked function and is periodic in Q with 
period b. For M large enough (see Fig. [4}, F(Q) is sharply peaked 
at Q = nb (n = 0, ±1, ±2, . . .) and with peak height of vM and 
width b/M. Thus for large M a useful approximation is to set F(Q — 
q) ~ ■\/~MdQ l q+nb- From Eq. l[4j we find that the spatial scale for 
this periodicity after expansion for time t is h as given in Eq. l[5). 
The validity conditions for d23b is the same as for the Wannier states 
discussed in the previous subsection. 

For later convenience we will give several useful relations for 
these basis functions that can be easily derived, 



]TF(Q-g)F(Q' T q) 

i 

F(Q - q)F(Q - q) 



MF(Q =f Q') 



F{Q-qf5 q 



(25) 



(26) 



The second relation expresses the delta-function like property of 
the quasi-momentum basis functions, however for small lattices (i.e. 
small M) this approximation may not be appropriate. 



III. CORRELATION FUNCTIONS 



Table I: Time scales and length scales for the experimental parame- 
ters for 87 Rb with V = 30E R , a = 425 nm, and M w 80. 



Xq 


Ma 


W(t) 


x M (t) 


t 


Ti 


T M 


58 nm 


34 am 


277 pm 


4.5 mm 


22 ms 


0.12 ms 


0.80 s 



A. Quantum field operators and correlation functions 

The many-body state of the system is described by the Hamilto- 
nian {j} in terms of the bosonic quantum field operator i/»( x ) which 
obeys the usual equal time commutation relations 



[^(x), ^i(x')] = 0. (27) 
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Correlation functions of the system can be expressed as expectations 
of products of the field operators at various positions and times (e.g. 
see f39ll4fll4lll ). Here we will mainly restrict ourselves to discussing 
normally ordered correlation functions at a fixed time (i.e. after ex- 
pansion). The correlation functions we will be most interested in are 
the expanded density: 



n(x) = <^(x),/>(x)}, 



(28) 



the first order correlation correlation function (one-body density ma- 
trix): 

G (1) (x,x') = (^(x)^(x')), (29) 
and the second order correlation function: 



G (2) (x,x') = ^ t (x)^ t (x')^(x')Vi(x)), 



(30) 



(see Ref. 14011 ). Normalized versions of these correlation functions 
can be defined as 



S «(x,x') = 
3 (2) (x,x') = 



GW(x,x') 
N /n(x)n(x') 
G (2) (x,x') 



(31) 



(32) 



n(x)n(x') 

Finally, we mention the density covariance function, defined as 

C(x,x') = (ft(x)n(x')) - (n(x))(n(x')>, (33) 

where the density operator is n(x) = i/^(x)?/>(x). This correlation 
function directly relates to the measurements made in experiments 
where density images are taken of the system and the shot noise an- 
alyzed. Commonly the normalized form of this observable is ana- 
lyzed, defined as 



g(x!x0= (n(x)n(xO) _ 1 
n(x)n(x ) 



(34) 



To allow us to evaluate these correlation functions after expansion we 
now turn to developing an expression for the field operator in terms 
of the single particle modes discussed in the previous section. 



B. Bose-Hubbard Hamiltonian 

In the lattice the ultra-cold system of bosons is well-described by 
the tight-binding Bose-Hubbard model 

<U> j j 

where J and U are respectively the tunneling matrix element be- 
tween sites and the onsite interaction |4"3l . and £j is the local potential 
offset of each lattice site (i.e. ej = Vcxt(Rj)). The operator a. cre- 
ates a boson in a ground vibrational state at site j = { j x , j y , j z } (i.e. 
at location Rj) of the lattice, and (i, j) in the summation indicates 
that the sum is restricted to nearest neighbours. 

Eq. ([35) can be derived from full Hamiltonian l[T) by considering 
the field operator projected to the ground band, i.e. 



(36) 



justified by the assumption that temperature and interaction effects 
are sufficiently small to remove the need to include modes of higher 
bands (e.g. see Refs. I42ll43ll ). In what follows we assume that prior 



to expansion the system is in a many-body eigenstate of Eq. i35\ . 
Of particular interest is that this ground state undergoes a quantum 
phase transition as the values of U and J change, for instance, by 
changing the lattice depth. For the ratio U/J below a critical value, 
g c , the state of the system is superfluid with long range phase co- 
herence. For U / J > g c the system enters the Mott insulating phase 
with a gapped excitation spectrum, suppressed number fluctuations 
and no long range phase coherence. In general g c depends on the 
dimensionality and the mean number of atoms per site. 



C. Field operators in the far-field limit 

We now consider the many-body expansion of the system when 
the external potentials are released at t — 0. In the Heisenberg pic- 
ture, the free evolution of the ground band field operator is given by 



(37) 



where we have taken the initial state to be entirely contained in the 
ground band (c.f Eq. d36t). Due to our approximation that the system 
freely evolves the time-dependence is completely contained in the 
single particle evolution (see Sec. Ill Bb . 

Under the condition that initially only the ground band is occu- 
pied, our expression for the field operator projected into the ground 
band l |37| l provides a useful and accurate description. However, care 
must be taken when evaluating correlation functions, as the commu- 
tation relation for ij) g (x) is non-local, i.e. 

[4(x,t),^(x',t)] = ^wo(x-Rj,i)«;o(x'-Rj,t). (38) 



For this reason it is more useful to work with normally ordered cor- 
relation functions, i.e. where annihilation operators precede creation 
operators. Indeed, n(x), G^(x,x') and G' 2 '(x, x') take the same 
form (i.e. Eqs. l|28t-l|30|l) whether we use the full or ground band pro- 
jected fields (assuming only ground band modes are occupied in the 
manybody state). On the other hand, the expression for the density 
covariance function in terms of ipg takes a different form compared 
to Eq. i33\ due to the non-local commutation relation (see Appendix 
lAl . However, it is more easily evaluated as 



C(x,x') = G (2, (x,x') +n(x)5(x-x' 



(39) 



In the above expression the normally ordered G* 2 ' (x, x') and n(x) 
functions can be evaluated using either ip or ip g . In the remainder 
of the paper we will concentrate on evaluating n(x), G^(x,x') 
and G (2) (x,x') using the ground band projected operators. These 
results can be directly related to the density covariance using Eq. l |39t 
if needed. 

Assuming that the expansion is sufficiently long for the far-field 
approximation to hold, we can express the expanded field operator in 
the simplified forms 

4(x,i) = J]a j A(x, i ) e ^ Q(x)R J, (40) 
j 

= £)& q j4(x,t)J*(Q(x)-q) J (41) 
q 

i.e. in Wannier and Bloch mode expansions respectively, 
where A(x,t) = A(x,t)A(y,t)A(z,t) and F(Q) = 
F(Qx)F(Qy)F(Q z ) are the three-dimensional generaliza- 
tions of the functions A(x,t) and F(Q) introduced in Sec. [II] The 
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operators 6 q are defined as the Fourier transform of the cij operators. 
We have suppressed time arguments on the mode operators which 
will be understood from hereon to take their t = (in situ) value 
and remind the reader that Q(x) also has an implicit dependence on 
time as indicated in Eq. {4j. 



IV. TRANSLATIONALLY INVARIANT RESULTS (ID) 

In this section we give the two-point first and second order corre- 
lation functions for a system of bosons in a translationally invariant 
(ej = 0) one-dimensional (j — > j) optical lattice of M = 15 sites 
for parameters corresponding to the superfluid and Mott-insulator 
phases. For consistency we will use the two sets of lattice param- 
eters given in Table Ull to realizes these phases, where the expansion 
time, t, is chosen sufficiently long for the far-field approximation to 
hold. For ease of comparison we will always plot our correlation 
functions against Q rather than x, since in terms of this variable of 
the correlation functions are independent of expansion time in the 
far-field limit. 



Table II: Time scales and length scales for the calculations presented 
in this paper. We take the case of 87 Rb, with a = 425 nm, and 
M = 15. 



0.1 





V 


t 


T M 


Xo 


Ma 


W(t) 


x M (t) 


SF phase 


9E R 


50 ms 


28 ms 


78 nm 


6.4 fim 


465 fim 


68 mm 


MI phase 


18 E R 


50 ms 


28 ms 


66 nm 


6.4 fim 


554 /im 


96 mm 



For compactness of equations in what follows we abbreviate our 
notation according to: 



A(x,t) 
A(x',t) 
F(Q(x) - q) 
F(Q(x') - q) 



A. 

A'. 

F Q-q, 



(42) 
(43) 
(44) 
(45) 



and so on. 



A. Pure Bose-Einstein condensate 

The first approach we consider is the prototype superfluid: A pure 
Bose-Einstein condensate. This approach is valid in the limit of van- 
ishing interactions (U = 0), but can also be used as an approximation 
for the interacting system in regimes where Gross-Pitaevskii theory 
is valid and quantum depletion can be neglected (see Sec. lIVDt . The 
ground state consists of all atoms in the 0-quasimomentum mode, i.e. 



3> 



(bl) N \0), 



(46) 



where TV is the total number of atoms. To calculate the correla- 
tion functions we substitute the quasimomentum expansion d41t of 
the field operator into Eqs. d28IH30|> for the correlations functions. 
Using that the only non-zero matrix elements are (&q&o) = TV and 
(SJSJSoSo) = N(N - 1), we obtain 



N\A["F, 



2 



n(x) — 
G {1 \x,x') = NA*A'F Q F Q ,, 
G {2 \x,x') = N{N -1)\A\ 2 \A'\ 2 F 2 Q F%, 



(47) 
(48) 
(49) 
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Figure 5: Comparison of expanded density distributions, (a) Perfect 
Mott-insulator (dotted line) and pure Bose-Einstein condensate (solid 
line) results, (b) Various superfluid phase results: pure Bose-Einstein 
condensate (solid line), Bogoliubov result (dashed line) and mean- 
field decoupling result (dash-dot line), (c) Various Mott-insulator 
phase results: perfect Mott-insulator result (dotted line) and includ- 
ing particle-hole corrections (dash-dot line). Parameters: TV ~ 85 
atoms in a lattice with M — 15. The Bogoliubov and meanfield 
decoupling calculations are for U /2 J ~ 3.2 and No ~ 65 with non- 
condensate atom numbers being TV ~ 20 (Bogoliubov) and N w 17 
(decoupling). The perfect Mott-insulator and particle-hole results are 
for U/2J w 42 with N « 85 (see Table[IIJl. 



These expressions depend on the number of sites M through the 
function Fq. In Fig. \5l a), and Figs.|6ja) and (c) we show examples 
of these correlation functions for the pure Bose-Einstein condensate. 

Result {47} is somewhat similar to the textbook example of laser 
light passing through a diffraction grating. Long-range coherence 
is revealed through the constructive interference that gives rise to the 
density peaks (see Fig.[5ja)). These peaks occur where Q is equal to a 
reciprocal lattice vector, as at these locations the phase accumulated 
by atoms propagating from distinct lattice sites are equal to within 
modulo 27T. In terms of the coordinate x the spacing between peaks 
is lb as given in Eq. l[5}. 

The superfluid correlation functions are (approximately) separable 
in the sense that 



G^(x,x') = \/ n(x)n(x'), 
G 2 {x,x) w n(x)n(x'), 



(50) 
(51) 



(to order 1/TV), with the consequence that the normalized correlation 
functions are 



g (x,x) 



A* A' 



g (2) (x,x') = 1- 1/TV. 



(52) 
(53) 
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Figure 6: Two point first order correlation function, G^ 1 ', for (a) pure 
Bose-Einstein condensate and (b) perfect Mott-insulator. Two point 
second order correlation function, G' , for (c) pure Bose-Einstein 
condensate and (d) perfect Mott-insulator. Parameters are the same 
as for the insulating cases considered in Fig. [5] 



separable. Indeed the correlation functions have well defined ridges 
running down the diagonal (Q = Q') and repeating periodically par- 
allel to the diagonals (i.e. Q = Q' +jb with j = 0± 1, ±2, . . .), see 
Figs-HJb) and (d) . The normalized first order correlation function is 



(x,x) 



A* A' F, 



Q-Q' 



\A\\A> 



iM 



(60) 



and reveals the absence of phase coherence with off-diagonal decay 
over a distance scale of Ax ~ ht/mM. 

The second order correlations have a diagonal ridge similar to 
those in G' 1 ', but also a smooth background envelope. The nor- 
malized form is 



(2)/ ' 

J (x,x 



1 + 



,(1)/ 



(61) 



where we have neglected terms of order n/N as being small. The 
second term (i.e. \g^'{x, x')\ 2 '), reveals the usual tendency for 
bosonic particles to cluster together. 

The results of this and the previous subsection emphasize the 
marked difference between the expanded correlations of the super- 
fluid and Mott-insulator phases in an optical lattice. 



C. Finite temperature ideal gas 



B. Perfect Mott-insulator 

In the limit where U 3> J, we can neglect the tunneling and ap- 
proximate the ground state of the Bose Hubbard Hamiltonian l |35| l as 
a perfect Mott-insulator 



M-l 



(54) 



3=0 



where n = N/M is the filling factor or mean number of particles per 
site. We have assumed that the number of particles is commensurate 
with lattice sites, i.e. that n is an integer. Using the Wannier expan- 
sion d41t of the field operator to evaluate Eqs. d28M30b we obtain 



n M (x) = N\A\ 



N 



A*A'F, 



M 

n\a\ 2 \a'\ 2 



Q- 



where we have used that 



{a\a]a k ai) 



(55) 
(56) 

,(57) 



(58) 



n (SikSji + SuSjk) — n(n + l)Sij5ki5jk, (59) 



and X) . e l( - Q - Q ' )R i = y/MF Q - Q ,. We have labelled these corre- 
lations with a subscript M for future reference. In Figs. |5]a), and 
Figs.(6{b) and (d) we show examples of these correlation functions 
for the perfect Mott-insulator. 

The expanded density of the Mott-insulator reveals the lack of 
phase coherence, i.e. the density arises as the incoherent sum of 
the individual Wannier states expanding from each site so that no 
interference peaks are observed (see Fig. 0a)). 

The correlation functions reveal many interesting features. Unlike 
the superfluid case, the Mott-insulator correlation functions are not 



It is of interest to compare the Mott-insulator properties with those 
of an ideal gas at finite temperature |44| , |4^ . |4^ . |47|1 . For this system 
the quasimomentum modes form a diagonal basis, which fluctuate 
according to 



{b ql b q 2) = n q i8 q i, q 2, 
(frqi&^VM = n q in q2 (l + 8 qhq2 ), 



(62) 
(63) 



in the grand canonical description, where n q is the mean occupation 
given by the Bose-Einstein distribution. The correlation functions 
are given by 



n(x) 
G (1) (x,x') 



G (2) (x,x' 



\A\ 2 J2n q F^ q , 
i 

A*A'Y^n q F Q ^ q F Q ,_ q , 
m 2 |Af Ew [Fq-^q'- 

qq> 

+F Q - q F Q i - q F Q , _ q , 7<Q_ 9 '] , 



\A\ 2 \A' 



(64) 
(65) 
(66) 

(67) 



2 

Q'-q 



In this treatment we have assumed that the system is not condensed, 
a point we address further in Sec. IV C| 



1. High temperature limit 

We now evaluate the correlation functions l!64b-(l67t in the high 
temperature limit, whereby the all quasimomentum states are occu- 
pied, however we assume that no higher bands are occupied (also see 
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14511 ). This limit could be arranged by taking a sufficiently deep lat- 
tice that £bw <C kT <C e g ap, where ebw ~ 4 J is the ground band 
width, and e gap ~ ft^Latt is the excitation gap to higher vibrational 
bands. In this regime the mean occupation of each of the quasimo- 
mentum levels is the same and is given by n q m N/M, i.e spatial 
filling factor n. The correlation functions now simplify to 



n(x) 
G w (x,x') 

G™(x,x') 



N\A\ 
N 



A*A'F Q _ Q , 



VM 
N\A\ 2 \A'\ 2 



N 



(68) 
(69) 

(70) 



where we have used l |25l > in deriving these results. Comparing these 
results with Eqs. d55t-(l57l> we see that the measured correlations of 
the high temperature ideal gas and the Mott-insulator state differ by 
corrections that are negligible when N is large. 

Aspects of the results presented so far have also been considered 
in Refs. fivl . r37L f38h . It is of interest to move beyond the pure Bose- 
Einstein condensate and perfect Mott-insulator approximations by 
considering: (i) the role of interactions in the superfluid state and 
(ii) tunneling in the Mott-insulator regime. In the following three 
subsections we turn to developing more realistic descriptions of the 
optical lattice system that go beyond those simple models and should 
be widely applicable to experimental regimes. 



D. Bogoliubov treatment of the superfluid phase 

The Bogoliubov treatment of the superfluid limit of the Bose Hub- 
bard Hamiltonian was given in Refs. Hmmnsni. 

and we 

refer the reader to those references for a detailed description of the 
method. 

Briefly, we recognize that condensation will occur in the 0- 
quasimomentum mode and set 



b ,i>l 



n%, No > i, 



(71) 



where No is the condensate occupation, and expand the Bose Hub- 
bard Hamiltonian to quadratic order in the remaining quasimomen- 
tum operators {6 9 ^o}- Performing the canonical transformation 



b q = u q (3 q - v q 0l , (g/0), 



(72) 



the Bogoliubov quasiparticle operators f3 q are introduced, where the 
amplitudes {u q ,v q } are taken to be 



e°+n U + ef 



2ef 



i rr B 

e„ + noU — e„ 



2ef 



(73) 
(74) 



with no = No/M the average number of condensate atoms per site 
and 



e q = 4Jsin 2 (ga/2), 



e°[e° + 2n C/]. 



(75) 
(76) 



This choice of transformation brings the quadratic Hamiltonian to 
diagonal form, i.e. H « const. + e f PqPq- 

Within the Bogoliubov approximation, the T — state of the sys- 
tem is given by the quasiparticle vacuum state, i.e. /3q^o|QP V ac) = 
0, and the condensate population is reduced from the total number 
of atoms due to the quantum depletion N, i.e. Nq — N — N with 



n = y 



J . These depleted atoms have an effect on the corre- 



lation functions as we now show. 

To calculate the correlation functions we transform the quasi- 
momentum expansion of the field operators J4 1 b to the quasiparticle 
form basis. Care needs to be taken with the four-field correlation 
function for which we use a number conserving approximation, i.e. 
the replacement fajfeo — N — X 9 ^o "qbq, rather than simply setting 
bjfeo = No (&!]. Following this procedure we obtain 



n(x) = \A\ 2 [N F^+J2v 2 q F ( 



Q-q] ' 



q^O 



G w (x,x') = A* A [NoFqFqi + ^2 v q F Q-i F Q 



q^O 



(77) 
(78) 



G< a >(*y) = \A\ a \Ar\[^+2^uivi]^,+^vUN -2ui)(F%F%,_ 9 +^_ q F' ( 



q^O 



\ 2 2 t-,2 r-,2 , \ ^ 4 P 2 J . \ ^ 2 2 j-,2 r? 

+ V q V q' F Q-q F Q'-q' + Z^ V 1 F Q-1 F Q'-1 + 1^ U 1 V 1 F Q-1 F ' 

q^O q^O 



,2 

Q'+q 



(79) 



In Figs. Ob) and|7ta) we show results for n(x) and G' 2 ' (a;, a;') in the Bogoliubov approximation. 
For this case it is interesting to consider the covariance which, according to Eq. l |39t , is 



C(x,x') = \Af\A\ 2 \ 2Y J ulvtFlFl,+Y J [vt F Q-q F l>-q+<4 F Q-q F l>+q] ~ 2 E U V* { F Q F Q'-q + F Q-q F ( 



k^O 



q^O 



q^O 



-\A\ 3 (N F^+^2v^_ q )S(x-x') 

q^O 



(80) 
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Figure 7: Second order correlation function, G^ 2 \ in the super- 
fluid regime, (a) Bogoliubov result and (b) meanfield decou- 
pling result, (c) Decoupling prediction for the pairing strength, 
G (2) (/f,/2, ~l h /2) (see text) as U/J is varied. The result for the 
case shown in (b) is indicated by a square. Parameters: Cases (a) and 
(b) are the same parameters as used in the superfluid calculations in 
Fig.0 



We note three observable effects of the quantum depletion in 
Eq. ( 180b . (i) The u^-term represents the additional fluctuations 
due to the quantum depletion which is of the same form as the 
thermal depletion for the finite temperature ideal gas. (ii) The 



u t v q^Q-q^Q' + 9 " term represents the pairing between parti- 
cles in the quantum depletion which, due to momentum conserva- 
tions, appears as a correlations between density fluctuations at +q 
and -q. (iii) The ~2J2 q ^u 2 q v 2 (F 2 F 2 ,_ q +F 2 _ q F 2 ,)-term 
represents the anti-correlation between depleted atoms and atoms 
in the condensate. This last term is absent if we do not make the 
number-conserving approximation. 

In Fig. [3a) the depletion and pairing effects can been seen. The 
weak diagonal line (e.g. along Q — Q') arises from the quantum 
depletion, whereas the other diagonal (e.g. Q — —Q') arises from 
the pairing in the system. This figure shows that the pairing and 
depletion are approximately the same size, and are much smaller 
than the signal due to the condensate which dominates at Q = or 
Q' = 0. Normalizing the correlation function to obtain g^ 2 '{x, x') 
significantly enhances the relative size of the pairing and depletion. 



E. Decoupling approach 

The meanfield decoupling approach is widely used to describe the 
many-body state of the system in both the Mott-insulator and the 
superfluid regions. The essence of this approach l52ll is the approxi- 

OjOk, with otj = (%), This turns 



mation a] fife 




-0.5 0.5 
Q [units b] 



-0.5 0.5 
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Figure 8: Second order correlation functions in the Mott-insulator 
regime, (a) Perfect Mott-insulator state and (b) particle-hole correc- 
tions to the perfect Mott-insulator state. The same parameters are 
used as for the Mott-insulator calculations in Fig. [5] 



the grand canonical form of Bose Hubbard Hamiltonian l |35t into the 
site-decoupled form f70Tl . i.e. 



(81) 



where N = Ylj o-j^j is the number operator, /i is the chemical po- 
tential [711 and 



Hi 



-a^ajaj 



-i + a j+ i) +h.c\ (82) 



+Ja*(ctj-i + Oij+\) 



The manybody ground state of d81b can be written in the (decoupled) 
product form 



3=0 



(83) 



where the is the ground state of Hj, and can be convenient ex- 
pressed in terms of the number states of a^aj. In the translationally 
invariant system we consider here, each site is equivalent and we can 
drop the j-label (i.e. setting ctj — > a) and solve the single site prob- 
lem Hj . 

The advantage of this approach is that it is able to represent num- 
ber states and coherent states, and thus contains as limits the pure 
Bose-Einstein condensate and perfect Mott-insulator states. How- 
ever, the decoupling approximation neglects correlations between 
sites beyond those described by the coherent coupling a, and thus 
can only provide a rather limited description of spatial correlations 
in the system. 

Within the decoupling approximation we obtain the following ex- 
pressions for the correlation functions 



n{x) = \A\ 2 [{N-M\ a \ 2 ) + M\a\ 2 F 2 Q ] 



G {1 \x,x) = A*A' 



N X c 



M n 



F Q _ Q , + M\a\ 2 F Q F Q , 



G (2 \x,x') = \A\ 2 \A'\ 2 [M X a+ M 2 \o.\*FqFq 

2 rr>2 



(84) 
(85) 
(86) 



+M{ X dFq +q , + XqFq-q') + M X b(F% + Fqi) 
+2AH F Q F Q ,(\a\ 2 xcF Q _ Q , + ?R{a 2 xh}F Q+Q/ )], 
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1. Comparison of decoupling and Bogoliubov approaches 



1 { a 




U/2J 

Figure 9: Behavior of meanfield decoupling parameters as lattice 
depth changes, (a) Commensurately filled lattice with n — 1 and 
(b) incommensurate case with n = 1.2. Results for a (solid line), j3 
(dashed line), 7 (dash-dot line) and 5 (dotted line) are shown. Mean- 
field estimate of the Mott-insulator transition for n = 1 is indicated 
with the vertical dotted line. 



where p = {hjOLj), 7 = (aja^a,), 5 = (a-a-ajCLj), are parame- 
ters that can be determined from the meanfield solution (i.e. \ipj) in 
Eq. (83)) and 



XA 


= M X 2 c+5- 


\0\ 2 + 2\a\ 2 (4n-3\a\ 2 ) 








-7*)}~2n 2 , 


(87) 


KB 


= (M-4)(| Q 


2 Xc-23?{q(q/3*-7*)}, 


(88) 


XG 


= n- \a\ 2 , 




(89) 


\D 


= P-ot , 




(90) 



where 5ft means the real part. In Fig. [9] we show the typical behav- 
ior of {a, P, 7, 6} parameters as the lattice depth increases. For the 
commensurate case [Fig. [9{ a )] we see that all of these parameters 
go to zero when the Mott-insulator transition occurs. I.e. when 
(U/J) > <?c, the local solutions, \ipj}, approach number states (if 
n is an integer). In the U/J — » limit, the local solutions 
approach coherent states with a — > \/n, P —* n, 7 — > 7i 3//2 and 
<5 — > n 2 . For the incommensurately filled lattice [Fig. [9fb)], the 
Mott-insulator transition does not occur and residual coherence ex- 
ists in the deep lattice limit. 

The above results for {a, P, 7, 5} allow us to understand how the 
parameters {xa, Xb,Xc, Xd} behave in the limiting cases: 

a. Superfluid limit: As U — > 0, we have 
Xa,Xb,XC,Xd — > and G' 1 ' and G^ agree with the ex- 
pressions of the pure Bose-Einstein condensate. (Note that in this 
limit the lower two lines of Eq. (86) vanish). In Figs. [5jb) and 
[Tib) we show decoupling results for n(x) and G^ 2 \x,x') in the 
superfluid limit. 

b. Mott-insulator limit: For ( U/ J) > g c (and n=integer) 
we have that xa — > Nn — n(n + 1) and xc —* n, while xb, Xd —* 
0. Under these conditions and G^ agree with the expressions 
of the perfect Mott state. (Again we note that in this limit the lower 
two lines of Eq. (86) vanish). 



It is of interest to compare the Bogoliubov and decoupling results 
in the superfluid regime to better understand the physics captured by 
each approximation (see Figs.|7Ja) and (b)). As the value of U/J 
increases both methods predict an increase in the quantum depletion, 
however for the Bogoliubov approach the depletion is colored ac- 
cording to Vg whereas for the Decoupling approach the depletion is 
uniform in quasi-momentum space and proportional to (n — |«| 2 ), 
i.e. xc- These differences are most apparent in the density predic- 
tions for each method shown in Fig. [5jb). 

While the Bogoliubov approach predicts no transition, and is 
hence invalid in the deep lattice limit, at a finite value of U/J the 
parameter a goes to zero and a purely incoherent interference pat- 
tern is produced. Similarly both methods predict pairing, colored 
and proportional to u 2 v 2 for Bogoliubov, while the decoupling re- 
sult is uniform and proportional to (P — a 2 ), i.e. xd- The uniform 
nature of the depletion and pairing predicted by the decoupling ap- 
proach indicates one of the main deficiencies of this theory: the de- 
coupling approximation makes all single particle excitations of the 
system (other than the q = condensate) degenerate. 



2. Strength of pairing correlations 

In view of the reasonable comparison of the decoupling and Bo- 
goliubov approaches, in Fig.[7jc) we show the strength of the pairing 
as predicted by the decoupling approximation. To do this plot the 
value of G (2) at the location x = l b /2 and x' =-4/2 (i.e. Q = b/2 
and Q' — —b/2). As the lattice depth increases the value of the pair- 
ing is observed to increase, saturating to a maximum value after the 
system passes through the Mott-insulator transition. 



F. Particle-hole corrections to the Mott-insulator phase 

Recent experiments done in 3D and 2D systems fl5ll25h have ob- 
served that phase coherence on short length scales persists even deep 
in the insulating phase. This behavior can be attributed to a coherent 
admixture of particle-hole pairs to the perfect Mott state for small 
but finite tunneling and can be estimated by using first order pertur- 
bation theory. We consider first a homogeneous system with filling 
factor n. In the limit of infinitely strong repulsion, U/J—> 00, the 
ground state is what we called a perfect Mott-insulator, Eq. (54)- As 
shown in Fig. |5]c) the expanded density distribution of this state is 
flat and consequently has zero fringe visibility. By treating the tun- 
neling term as a perturbation to the interaction term, one can account 
for finite tunneling corrections and to first order in J/U one gets 
that on top of the Mott-insulator core the system has an admixture of 
nearest-neighbor "particle-hole" excitations (a site with an additional 
particle while one of its nearest neighbors is missing a particle): 



|$Ml}«|$M>+yXX ta '-|* M >> (91) 
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The restoration of the short range phase coherence induced by such 
excitations is signaled in the correlation functions which become: 



"Ml — nyi 



A 2 2N(n + l) 
U 



E(Q), 



(92) 



U M1 



r<( 2 ) - 



G® - ME^±A Fq _ qi{E{Q) + E(Q'U93) 



r>(2) _ 



AMa 2MN[E(Q)+E(Q')Y ] . 



U 



i(n + l)(l + 



F Q-Q' ■ 
M ■ 



M 



(2n 2 + 3n + 1) 



where E(Q) — —2Jcos(Qa) is the single particle dispersion re- 
lation. In Fig. [Sjc) the density «mi (x) is shown while in Fig. [8ja) 
and (b) the second order correlation functions for the perfect Mott- 
insulator and the case including particle-holes corrections are com- 
pared. A rather significant effect of these particle-hole correlations 
are observed as long wavelength modulations of the density, charac- 
terized by the term E(Q) given above. 

While the particle-hole modulation [see Fig.[5]c)] to the expanded 
density distribution has been experimentally measured fi"5l [2^1 the 
corrections to the noise correlations has not yet been observed. One 
possible reason is that in the noise correlation experiments part of 
the superfluid component had to be masked out before computing 
the correlation functions. 

Even though here we will limit our analysis to first order cor- 
rections we conclude this subsection by briefly outlining the basic 
physics beyond it: Recent theoretical analysis based on the Random 
Phase Approximation (RPA) [530 supports a physical picture of the 
system as a (dilute) gas of partice-hole pairs, mobile through the lat- 
tice, on top of a regularly filled Mott-insulator. Consequently what 
happens for larger J/U values is the development of phase coher- 
ence at longer length scales as the particle-hole pair size extends over 
larger distance scales than a single lattice constant through higher- 
order tunneling events which become non-negligible. 



V. GAUSSIAN APPROACH (ID) 

In the previous section results for the expanded correlations of a 
one-dimensional system of bosons in a translationally invariant lat- 
tice were derived. Extending these results, particularly for G^ 2 \ to 
higher dimensions and inhomogeneous potentials is difficult. Here 
we develop a different approach: we approximate the expanded mat- 
ter wave field as being Gaussian, so that all higher order correla- 
tion functions can be determined from the low order moments, i.e. in 
terms of the condensate wavefunction, and the normal and anoma- 
lous density matrices (introduced below). This approach avoids di- 
rectly evaluating the four-point correlation function and is thus better 
suited to numerical calculations, which are generally required for in- 
homogeneous situations (which we consider in the next section). We 
emphasize that while a Gaussian approach is in general insufficient 
for capturing the in situ many-body physics of the lattice system, 
particularly in the Mott-insulator regime, we show that it provides 
a good representation of the low order noise correlations of the ex- 
panded system for the various approaches we have considered so far. 

The general Gaussian expression for the G^ 2 \x,x') correlation 
function is 

G i2) (x,x') « \G (1) (x,x')\ 2 +n(x)n(x') - |$o(z)| 2 |$o(a:')| 2 
+ {$* {x)$* {x')M{x,x')+c.c} (95) 



where we have introduced the following quantities: 

Condensate orbital: $o (x). When the system is in the superfluid 
phase this quantity is non-zero and reflects the long-range order in the 
system. The condensate orbital can be obtained as an expectation of 
the field operator <&o{x) = (ip g (x)) in symmetry broken approaches 
or by the Penrose-Onsager criterion l54ll . 

Anomalous average: M(x, x') (or anomalous density matrix). In 
symmetry broken approaches, when a condensate exists in the system 
a fluctuation operator can be defined as 

£(x)=i> a (x)-Mx), ( 96 ) 
such that (£(x)) = 0. The anomalous average is then given as 



(97) 



The first line of Eq. ^ was used in Ref. |40|| to analyze the in 
situ correlations of a harmonically trapped Bose gas. Here we use the 
more general expression which also includes anomalous terms (e.g. 
see Chapter 4 of Ref. l55ll ). We reiterate that here we are not approx- 
imating the in situ Bose-Hubbard state as being Gaussian, only the 
matter wave field after expansion. Indeed, of the various approaches 
we have considered, only the pure superfluid, ideal gas and Bogoli- 
ubov methods are themselves Gaussian descriptions of the many- 
body state. The perfect Mott insulator, decoupling and particle-hole 
approaches are non-gaussian states and for these cases it is not clear 
that expressing the G^ function in terms of lower order moments 
l |95t is appropriate - an issue we examine more closely in the remain- 
der of this section. 

More generally, the expanded correlation functions can be ex- 
pressed in terms of cumulants which, as opposed to operator mo- 
ments, tend to become smaller at high orders |[56ll . Indeed, for a 
system described by a Gaussian density matrix, all cumulants higher 
than second order (which we refer to here as first order correlation 
functions) are zero j72j[. This approach would allow one to systemat- 
ically go beyond the Gaussian approach by considering higher order 
cumulants, however our results here suggest this is unnecessary. We 
also note that a more precise number conserving definition of the 
anomalous average could be adopted l57ll . however for the purposes 
of this paper the definition given in l |97t is adequate. 

We now apply the Gaussian approach to the ID translationally 
invariant cases considered in the last section to justify its use. 



A. Gaussian pure Bose-Einstein condensate 

For the pure Bose-Einstein condensate state considered in 
Sec. |IVA] we find that $ (x) = AVnFq and M(x,x') = 0. So 
that from Eq. ( 1951 ) we obtain 



G (2) {x,x') = N 2 \, 



\ 2 \A'fF%F%, 



(98) 



which differs from our earlier result J49t only by the N(N — 1) co- 
efficient on G^ which is negligible when N is large. More gener- 
ally, in this limit the system is completely described by the coherent 
condensate mode and all higher order correlation functions trivially 
factorize. 



B. Gaussian perfect Mott-insulator 



+|M(a 



For the perfect Mott-insulator state considered in Sec. IIVB1 the 
condensate and anomalous average are both zero and from Eqs. d55t . 
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J56b and {95j we obtain 

G [2) (x,x) = N\A\ 2 \A' 



N 2 

N+ M F ^' 



(99) 



This differs from the exact result l !57t by a correction to the first N 
in the square brackets that is negligible for large N. 



C. Gaussian finite temperature ideal gas 



For the finite temperature ideal gas considered in Sec. HVCl the 
Gaussian approach gives identical results to those in Eq. {67) [hence 



we also obtain the same predictions for the high temperature limit 
given in Eq. J70H. We note that in obtaining these results we assumed 
the condensate orbital and the anomalous average are zero. The no- 
condensate assumption was also assumed in the derivation presented 
in Sec. lIVCl whereby all the modes were taken to fluctuate according 
to the grand canonical ensemble. When a condensate is present the 
&o(x) corrections in Eq. l |95t are important (this issue is discussed 
more fully in Ref. |40|] ). In this case the q = mode defines the 
condensate, i.e. §o(x) = A-^/WoFq, and in the ideal limit (where 
there is no anomalous correlation) we obtain 



G {2) (x,x') = \A\ 2 \A'\ 2 \ 5> 9 rv [fQ- 9 fQW +F Q . q F Q ,_ q F Q ,_ q ,F Q _ q ,} -n 2 F%F ( 



Q' ( ' 



9 l' 



Y n i n i' F Q^i F Q'-i' + ^ n q FQ_ q FQ,_ q 

q q' q^O 



(100) 



(101) 



D. Gaussian Bogoliubov treatment of the superfluid phase 



The Bogoliubov treatment given in Sec. lIVDl is essentially a Gaussian theory and is well suited to the development of this section. The 
condensate orbital is <&o{x) — Ay/ TJqFq, so that the fluctuation operator is given by £(x) = ^E,/oK^ — v qfl- q )FQ- q and the 
anomalous average is 

M(x, x) = —AA' u q v q F Q ^ q F QI+q . (102) 

qjtO 

Thus using these results and those in Eqs. J77K l |78t and l |95t , we obtain 



G {2) {x,x') = 



N 2 FqFqi + Y v 2 q N {FlFl^ q + Fq_ 9 Fq 



q^O 



+ J2 V q V q> F Q-q F Q>-q> + J2 V t F Q ~ 1 F Q' ~1 + J2 U ^ F Q ~1 F Q' +1 
qq'^0 9/0 q^0 



(103) 



Compared to Eq. l |79t this result lacks the anti-correlation between the condensate and quasi-particle modes, but otherwise shows the same 
behavior for the quantum depletion and pairing. The pairing term arises from the anomalous average, and shows the importance of including 
the anomalous average in any description correlations for the interacting superfluid regime. The anti-correlation effects could be recovered by 
using a number-conserving generalization to the definition of the anomalous average, but we will not address this further here. 



E. Gaussian form of particle-hole corrections to the Mott-insulator phase 



Using the Gaussian approximation, one can calculate the particle-hole corrections to the perfect Mott-Insulator state considered in Sec. lIVFl 
by using Eq. l |95l > with the results given in Eqs. l |92| > and d93t . and by setting both the condensate and anomalous average to zero. This procedure 
yields, to first order in J/U, the following expressions for the correlations functions: 



G^ 2 \x,x') = |G^(a;,a;')| 2 + n M (x)n M i(x') + 25R[G[J ) (x, x')G^\(x, x')] + n M (x)n M i(x') + n M i(x)n M i(x'), 



N\A\ 2 \A' 



N + ^F ( 



M' 



Q-Q' 



\A\ 2 \A 



, l2 2 MN{E(Q)+E{Q') 
U 



n(n + l) 1 + 



(104) 



(105) 



Similarly to the perfect Mott-insulator situation, the Gaussian approximation results differ from the exact results, Eq. 
square brackets of Eq. J105t that is negligible in the thermodynamic limit. 



M 

|, by a term in the last 
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F. Gaussian decoupling approach 

The mean-field decoupling approach was introduced in Sec llV hi The condensate orbital is given by &o(x) = (ip g (x)} = ^/MaAFQ, so 
that the fluctuation operator is given by = ip g (x) — &o(x) and 

M(x,x') = AA'VM(f3-a 2 )F Q+Q ,. (106) 
Using this result and Eqs. d84b . J85t and l |95t we obtain 

G^(x,x') = \A\ 2 \A'\ 2 [MxA+M 2 \a\ A F%F%,+M(xlF% +Q ,+ X cF%- Q >) 

+M X b{F% + Fqi) + 2M% FqFqi {\a\ 2 Xc F Q . Q , + $t{a 2 X * D }F Q+Q ,)] , (107) 

I 




2 4 6 8 10 12 14 

U/2J 



Figure 10: Comparison of Gaussian and full decoupling results, (a) 
Comparison of xa for n = 1 (stars) and n = 1.2 (crosses), and xa 
for n — 1 (circles) and n — 1.2 (dots), (b) Comparison of xb for 
n = 1 (stars) and n — 1.2 (crosses), and x_b for n = 1 (circles) and 
n = 1.2 (dots). Meanfield estimate of Mott-insulator transition for 
?i = 1 is indicated with the vertical dotted line. 

with xc and xd as defined in Eqs. l |89t and l |90t , and 

XA = M X c, (108) 
Xb = M\a\ 2 Xc - (109) 



the Mott-insulator transition occurs and beyond this we have xa — > 
n(N — n — 1) and xa — ► nAT which agree to order n/N. For the 
incommensurate case no sharp transition is observed, however xa 
and xa are seen to be in reasonable agreement for all values of the 
lattice depth. 

Second, we compare xb and x.b in Fig. UOi b). For the com- 
mensurate lattice both quantities approach the value of in the pure 
Bose-Einstein condensate and in the Mott-insulator limits. In be- 
tween those cases, in the interacting superfluid regime, xb and xb 
take non-zero values and are similarly behaved. For the incommen- 
surate case xb and xb have appreciable values in the deep lattice 
limit. 

These comparisons give us confidence that the Gaussian decou- 
pling approach provides a reasonably accurate representation of the 
full decoupling results. We also note that in general the agreement 
improves as the filling factor increases. 

More generally, the results of this section strongly support that 
a Gaussian representation of the expanded correlation functions is 
accurate even for non-Gaussian states many-body states, such as the 
perfect Mott insulator, decoupling state and particle-hole result. For 
example, in the decoupling approach the single site moments, a, /3, 7 
and S are all important and poorly approximated as Gaussians (i.e. 
approximating 7 and S in terms of products of a and /3). However, 
the far-field matter wave results from the interference of atoms from 
all parts of the original lattice. In this case the central limit theorem 
ensures the measured statistics (in momentum space) are closer to 
Gaussian than were those of the in situ state. This suggests that the 
general applicability and accuracy of the Gaussian approach should 
be quite broad if we can accurately calculate the low order moments. 



VI. INHOMOGENEOUS GAUSSIAN THEORY (3D) 



1. Comparison with full decoupling result 

Here we compare the differences between the Gaussian decou- 
pling result above with the full result developed in Sec. IIVEI As 
emphasized in our choice of notation above, these differences arise 
between xa [Eq. d87H and xb [Eq. d88H and the respective quanti- 
ties xa [Eq. j!08H and xb [Eq. j!09H . These quantities are com- 
pared for cases of commensurate and incommensurate filled lattices 
in Fig. [10] 

In Fig. UOf a) xa and xa are considered. In the pure Bose-Einstein 
condensate limit (i.e. small U/2J) both results agree with xa — > 
and xa —> 0. The deep lattice (i.e. large U /2 J) behavior depends on 
whether the system is commensurately filled or not. For n=integer 



It would be feasible to extend the exact evaluations of G*- 2 ' con- 
sidered in Sec. lIVI to the 3D case for a translationally invariant lattice. 
However, in the presence of additional inhomogeneous potentials (as 
in experiments), the generic exact results for the second order cor- 
relation function would involve computing the four- point correla- 
tion functions ^(ajajafcfi;)^ without the simplifications afforded by 

translational invariance. In the context of 3D calculations, such re- 
sults would be computationally impractical to deal with. 

The primary advantage of the gaussian approach we have justi- 
fied in the previous section is that it offers a convenient method for 
extending the various techniques for calculating correlations to in- 
homogeneous 3D situations, and only involves 2-point correlation 
functions. This ultimately lends itself to a numerically tractable the- 
ory for calculating results that can be compared with experiments. 
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In this section we consider a full 3D formulation, and thus extend 
our previous abbreviated notation for this purpose. In particular, in 
this section we will use the following: 



A(x,t) -> A, 
A{x!,t) - A', 
F(Q(x)-q) - F q _ q , 



(110) 



and so on. 

For application to the inhomogeneous system it is useful to define 
the momentum space representation of system quantities. In this sec- 
tion we will indicate Fourier transformed quantities with a tilde, e.g. 
the discrete Fourier transform of the density yields 



1 

— W / n i e 
Mi V 



q.Rj 



(111) 



where q is limited to the discrete set of allowed quasi-momentum 
values. 

It is also convenient to extend this definition to the continuous 
Fourier transform which we can parameterize with far-field momen- 
tum Q. We can easily show that 



riQ 



Q.Rj 



1 

Mi 



Q-qi 



(112) 



i.e. the Fourier transformed density (indicated with a subscript Q) 
contains the Fq-q implicitly, and thus using this form affords us 
a more compact notation. We emphasize that this quantity is not 
the momentum density, which could be obtained from the Fourier 
transform of the field operator. In many of the derivations that follow 
we also make use of the identity 



q-pQ'-q' 



(113) 



A. Inhomogeneous pure Bose-Einstein condensate 

The pure Bose-Einstein condensate in the inhomogeneous system 
is described by the discrete Gross-Pitaevskii equation 



fjbZi = —J^2 z i + e ' z i + t^|zi| 2 Zi, 

(ij> 



(114) 



where z\ = (ai) is the condensate amplitude on site i, and /i is the 
chemical potential chosen to ensure \z\\ 2 = N. Applying the 
Gaussian formalism requires the following quantities 



n(x) = |A| 2 Af 3 |2 Q | 2 , 
<E> (x) = AM^zq, 
G (1) (x,x') = A*A'M 3 zqz q , 



(115) 
(116) 
(117) 



with the anomalous average taken to be zero. We reminder the 
reader that zq is defined as the Fourier transform (see above) of 
zi. For the results in this section we do not give G' 2 '(x,x') since 
it is (within the Gaussian approximation) completely determined by 
{$o(x), n(x), G (1) (x, x'), M(x, x')} using {95}. 



I 
o 
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Figure 11: Comparison of G^ 2 ' correlation functions for ID lattice 
systems: (a) Superlattice with periodically modulated density, (b) 
uniformly filled translationally invariant lattice. Both cases calcu- 
lated using the inhomogeneous perfect Mott-insulator approach (i.e. 
tunneling is neglected) with N = 29 (superlattice), N = 30 (uni- 
form lattice) and M = 15. 



B. Inhomogeneous perfect Mott-insulator 



Neglecting tunneling, we take the manybody state to be of the 
product form 

i* M )=n_i= (<)">>, (us) 



where nj = [(m ~~ e j)/U\ is the integer site occupation, with /i the 
chemical potential determined by the condition that ftj = N. 

For this state the condensate fraction and anomalous average are 
zero and we obtain 



G 



n M (x) = |yl| 2 iV, 
(1) (x,x') = A*A'M%nQ-Qi 



(119) 
(120) 



This result shows that first and second order correlation functions 
reveal information about the Fourier transform of the density in the 
system. 



1. Application to a super-lattice system 

We now develop a simple application of the inhomogeneous per- 
fect Mott-insulator formalism. We consider a ID system in which an 
external potential (i.e. the tj) causes the many-body state to be of the 
form given in Eq. l l 1 I St with 



3, j even 
1, jodd 



(121) 



i.e. a system with a modulated filling. Such an density arrange- 
ment can be produced using a superlattice formed by superimposing 
several optical lattices (e.g. see l58l [59ll ). In Fig. II lf a) we show 
the second order correlation function for this state [compared to the 
uniformly filled Mott-insulator state shown in Fig. II l[ b)l. The addi- 
tional spatial structure due to the density modulation is clearly appar- 
ent through the introduction of a new periodicity scale at wavevector 
b/2, i.e. a length scale corresponding to the 2-lattice site modula- 
tion in the filling density. Similar conclusions have been reached in 
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a recent experiment analyzing the noise correlations in a two-color with nj = nj — |aj| 2 the incoherent density, and (o? 2 )q the Fourier 
lattice [161. 



transform of aj . 



C. Inhomogeneous ideal gas 

The ground band modes of the (ideal) inhomogeneous Bose Hub- 
bard Hamiltonian satisfy 



E s u\ = — J Mj + eiuf, 

<ij> 



(122) 



where E s is the energy eigenvalue and uf is the amplitude on the i-th 
site of the s-excited mode. For consistency with the other treatments 
in this section we take the s = mode to be that of the condensate 
(when present) and label its amplitude as y^ou? — > zi. We obtain 
that 



$o(x) = AM^zq, 
n(x) = \A\ 2 M S l\2 Q \ 2 +J2 



n s \ii a Q \ 2 



(123) 
(124) 



G (1) (x,x') = A*A'M 3 yz& Q , +Y^n s u%u s Q ,j ,(125) 

where n a is the mean thermal occupation of the a mode and the 
anomalous average is zero. 



E. Inhomogeneous Bogoliubov approach 

The extension of Bogoliubov theory to inhomogeneous systems 
has been given in Refs. l44l Boll . Briefly, the discrete Gross- 
Pitaevskii and Bogoliubov equations 

i + t/h| 2 ] 2j , (130) 
J ^ <4 + [Vj cfi - (i]uj - Uz?v-, (131) 





— J ^2 2m 

<mj> 




-j M ™ 
w 




-J 5^ "™ 

<m,j> 



need to be solved to describe the system, where Vj cff = ej + 2[/|zj | 2 
is the effective potential and s labels the quasi-particle excitations 
described by the modes {it? , v? }. 

Using these solutions we determine the condensate mode as 



$ (x) = AM 2 zq, 
and the fluctuation operator as 

£( x ) = AM* ^(flqft — Vqfil) 



(133) 



(134) 



We also obtain that 
n(x) 

G (1) (x,x') 
M(x,x') 



\A\ 2 M 3 (\z Q \ 2 +^\v Q A, 

A*A'M 3 ^z* q z q , +J2*Q*<* 
—AA'M 3 ^2 uqVq' • 



(135) 

(136) 
(137) 



D. Inhomogeneous meanfield decoupling 

The meanfield decoupling approach is immediately extensible to 
3D from the form presented in Sec. lIVEl bv introducing the appropri- 
ate generalizations to the parameters introduced earlier, i.e. rii — * ni, 
oti — > Qi, and /3i — > In this approach the condensate is identified 
as the mean value of the field operator, and we obtain that 

$o(x) = AM§&q, (126) 
n(x) = m 2 {Mi,^ + A/ 3 |a q | 2 }, (127) 

G (1) (x,x') = AM'{M*nJ,_ Q , +A/ 3 5q5 Q ,}, (128) 

M(x.x') = AA'{Mf% Q ,-Mi(?) Q+Q ,},(129) 

I 

Using perturbation theory to account for finite tunneling corrections to the inhomogeneous extension of the perfect Mott-insulator state 
(Eq. d 1 1 8b >. again just as in Sec. lVIBl the the condensate fraction and anomalous average are zero and we obtain 



F. Particle-hole corrections to an inhomogeneous Mott 
Insulator 



Jala; 



|*M1> « |*M> + V J-— " 

t-f U(l + m - r 

<»J) 



-njj+ei-ej 



I*m), 



(138) 



"Ml 



for the first order perturbed wavefunction (with |3>m) given in dl 18b ) and the following expressions for the correlation functions 

J(rii + l)rij 

J(n ; + l)nj 

<i.j> 



(x) = n M (x) + 21^1 2 ^cos (Q ■ [R| - R,]) — 



U ( 1 + m — rij ) + ei — ej 



(139) 
(140) 



(7(1 + 7li — rtj ) + £i — £j 

G«(x,x') = GiJ ) (x,x')+^'E e * (Q ~ Q ' ,Ri ' Q ' (i " • ' <J " Rjl 

I 

VII. CONCLUSIONS AND OUTLOOK released from an optical lattice. We have given validity conditions 



The central purpose of this paper has been to develop a formalism 
for describing the correlations observed between ultra-cold bosons 
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for a far-field regime, which allowed us to simply relate the corre- 
lations measured after expansion to the in situ correlations of the 
system in the lattice. Since the in situ state is well-described by the 
Bose-Hubbard model, we have developed our formalism for a va- 
riety of commonly used methods for solving this model, including: 
the Bogoliubov method, the meanfield decoupling approach, and the 
particle-hole perturbative solution about the perfect Mott-insulator 
state. These methods provide a rather general set of solutions appro- 
priate to a broad range of behavior in the lattice system. We have 
also considered a Gaussian version of our formalism that gives gen- 
eral correlation functions of the expanded system in terms of the first 
order correlation functions. We justify this approach by compari- 
son to our earlier (non-Gaussian) results and show that it provides a 
complete and accurate description of the results in all regimes consid- 
ered. Finally, we used the Gaussian approach to develop a tractable 
3D formalism appropriate to inhomogeneous systems - i.e. relevant 
to current experiments. The application of this formalism to a variety 
of experimental regimes will be the subject of future work. 

So far experiments have only analyzed the basic correlation prop- 
erties of bosons in the Mott-insulator regime including second order 
correlations EQSGJ] and density profiles fH[Tlr3l. Technical 
noise has so far made measurements of the superfiuid phase difficult, 
so it would be of interest to see continued refinement of experimen- 
tal techniques into this regime. Possibly the most interesting region, 
from the theoretical perspective, is the nature of correlations in the 
transition region. Preliminary studies of some correlation aspects 
near the transition have been considered in Ref. fl5ll . 

Our formalism now presents numerous avenues for future theo- 
retical development. A few such areas of interest are: (i) Char- 
acterizing correlation function behavior in the region of the super- 
fluid to Mott-insulator transition, (ii) Investigating the influence of 
temperature on interacting cases, both in the superfiuid and Mott- 
insulator regimes, (iii) Considering the influence of disorder on 
the system phase diagram (e.g. emergence of the Bose-glass phase 
03, HO, |6lL [62|, l63ll ) and how signatures of these phases might 
emerge in correlation functions. Some of these problems will re- 
quire approaches that go beyond those covered here (e.g. may re- 
quire Monte-Carlo solutions [64J,.65j]), however our results suggest a 
Gaussian approach will be suitable, i.e. only requiring knowledge of 
{$o(x),G^(x,x'),M(x,x')}. 
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Appendix A: PROJECTED DENSITY COVARIANCE 

Here we examine the properties of the density covariance func- 
tion constructed from the projected field operator. To be precise we 
decompose the full field operator into the form 

^(x,t) = ^ 9 (x,t)+4(x,*) ) (Al) 



where i/> 9 (x, t) is the ground band operator defined in Eq. d36t and 
tl> e (x,t) describes the modes orthogonal to the ground band. The 
commutation relations for these operators are 

|&(x,t),$(x',t)] = A(x,x',i), (A2) 

|&(x,t),$(x',t)] = 5(x-x')-A(x,x',t), (A3) 

[&(x,t),#(x',t)] = 0, (A4) 

where A(x, x', t) = J]V u>o(x — Rj , t)luo(x' — Rj, t). 

The experimental observable is the density covariance as defined 
in Eq. d33t . Expressed in terms of xp g and i/> e , we obtain 

C(x, x') = C a (x, x') + G (1) (x, x') {5(x - x') - A(x, x', t)}(A5) 

where we have introduced the ground band density covariance 

G 9 (x,x') = <n 9 (x)n 9 (x')) - (n 9 (x))(n 9 (x')>, (A6) 

with n 9 (x) = ^> 9 (x)i/S 9 (x), and have assumed the excited bands to 
be in their vacuum state so that (x')i/S e (x')) = 0. 

Appendix B: FINITE RESOLUTION 

The issue of finite imaging resolution and the effect of taking a 
column density has been dealt with in reference fI3l . however we 
review the basic considerations here for completeness. 

For the case of absorption imaging, the measured density is given 

by 

h m (x,y) « j dx' dy' PSF (x' ,y') J dzn(x — x',y — y',z), (Bl) 

where the integral along z accounts for taking column density and 
the point spread function, PSF(x,i/), accounts for the limited opti- 
cal and detector resolution. The point spread function can be ap- 
proximated as a normalized Gaussian function of widths a x and 
tj y along the x and y directions respectively, i.e. PSF(s,j/) w 
ITj^a: y {exp(—x 2 /2a 2 )/y/2naj}. We note that for the case where 
the atoms are counted directly, e.g. using a micro-channel detector 
plate in a meta-stable helium system, we can associate a 3D point 
spread function for the measurements which includes the (2D) spa- 
tial resolution of the detector and the time resolution. 

The general situation of interest is when a peaked feature occurs 
in a correlation function. To understand the effects of imperfect res- 
olution we consider a ID example to develop the basic ideas. Take 
the peaked function to be of the form 

/peak(a;) = he X p[-x 2 /2(Axf], (B2) 

with feature height h. After convolution with the (ID) point spread 
function PSF(x) = exp(— x 2 /2a 2 ) /y/2~na x we obtain the mea- 
sured quantity 

fpclt(x) = hA x exp (--^y^J , (B3) 

where 

A. = 1 (B4) 

W x = V(Aa;) 2 + al, (B5) 

are the resulting height ratio and width parameters respectively. 

Several approximate limits for the effects of imaging on the height 
and width of the features in f™ca£{ x ) are easily found: 
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a. perfect resolution: Ax 3> a x 

In this limit the resolution is much better than the characteristic 
size of the peak feature and we have 

W x w Ax, (B6) 
A, « 1. (B7) 

In this case the features of the peak are accurately revealed by the 
imaging. 



smallest feature size (Ax) m i n = h/M where If, = hbt/m is the 
periodic length scale of the expanded system (also see Eq. ®). In 
practice the largest value of the imaging resolution is (<7j) ma x = 
lb, since beyond this length scale the underlying periodicity in the 
system will cause the peak to reoccur. Thus, in the no resolution 
limit we have that (from Eq. l |B9t ) 

A, « -!-. (BIO) 



b. limited resolution: Ax <C o~ x 

In this limit the resolution is larger than the characteristic size of 
the peak feature and we have 

W x « a x , (B8) 
A, « — < 1. (B9) 

<Tx 

Thus the measured peak feature is strongly altered by imaging, in 
particular the width is that of the imaging resolution and the feature 
height is reduced according to the ratio of the feature size to imaging 
resolution. 



c. no resolution 

For the lattice system there is usually a limit to the ratio Ax / a x 
appearing in Eq. (B9\ . The finite number of lattice sites makes the 



1. Generic experimental situation 

Generally experiments (e.g. Refs jTJ 03, EUl ) have operated in 
the regime where Axj w h/M <C Oj <§C h for j = x, y (in the 
imaging plane) and taking a column density (i.e. no resolution limit) 
along 2. Most importantly this means the height of peaks in the cor- 
relation function are reduced by a factor of 



= — ■ 



A more complete treatment yields an additional factor of 1/47T f73h . 
The widths of the peaks will be given by o x and o y in the x and y 
directions respectively. 
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